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Plasma turbulence described by the Hasegawa-Wakatani equations has been simu- 
lated numerically for different models and values of the adiabaticity parameter C. It 
is found that for low values of C turbulence remains isotropic, zonal flows are not 
generated and there is no suppression of the meridional drift waves and of the par- 
ticle transport. For high values of 6, turbulence evolves toward highly anisotropic 
states with a dominant contribution of the zonal sector to the kinetic energy. This 
anisotropic flow leads to a decrease of a turbulence production in the meridional sector 
and limits the particle transport across the mean isopycnal surfaces. This behavior 
allows to consider the Hasegawa-Wakatani equations a minimal PDE model which 
contains the drift- wave/zonal- flow feedback loop prototypical of the LH transition in 
plasma devices. 
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I. INTRODUCTION 



One of the major experimental discoveries in nuclear fusion research was the observation 
of a low-to-high (LH) transition in the plasma confinement characteristics^. This transition 
results in significantly reduced losses of particles and energy from the bulk of the magnet- 
ically confined plasma and, therefore, improved conditions for nuclear fusion. Since this 
discovery, LH transitions have been routinely observed in a great number of modern toka- 
maks and stellarators, and the new designs like ITER rely on achieving H-mode operation 
in an essential way. The theoretical description of the LH transitions, and of nonlinear and 
turbulent states in fusion devices, is very challenging because of the great number of im- 
portant physical parameters and scales of motion involved, as well as a complex magnetic 
field geometry. Direct numerical simulations (DNS) of gyrokinetic Vlasov equations have 
become a popular tool for studying such fusion plasmas, which involves computing particle 
dynamics in a five- dimensional phase space (three space coordinates and two velocities) and, 
therefore, requires vast computing resources. Physical mechanisms to explain the LH tran- 
sition have been suggested. One of these mechanisms is that small-scale turbulence, excited 
by a primary {e.g. ion-temperature driven) instability, drives a sheared zonal flow (ZF) via 
a nonlinear mechanism, through an anisotropic inverse cascade or a modulational instabil- 
ity. After this, the ZF acts to suppress small-scale turbulence by shearing turbulent eddies 
or /and drift wave packets, thereby eliminating the cause of anomalously high transport and 
losses of plasma particles and energy. 

Importantly, such possible scenario to explain the LH mechanism was achieved not by 
considering complicated realistic models but by studying highly idealised and simplified 
models. More precisely, generation of ZF's by small-scale turbulence was predicted based on 
Charney-Hassegawa- Mima (CHM) equation^'^ very soon after this equation was introduced 
into plasma physics by Hassegawa and Mima in 1978^, and even earlier in the geophysical 
literature-'. The scenario of a feedback in which ZF's act onto small-scale turbulence via 
shearing and destroying of weak vortices was suggested by Biglari et al. in 1990^ using 
an even simpler model equation, which is essentially a 2D incompressible neutral fluid de- 
scription (equation (1) in Refj6|). Probably the first instances where the two processes were 
described together as a negative feedback loop, turbulence generating ZF, followed by ZF 
suppressing turbulence, were in the papers by Balk et al. 1990^'*^. Balk et al. considered 



the limit of weak wave-dominated drift turbulence, whereas the picture of Biglari et al. ap- 
plies to strong eddy-dominated turbulence. In real situations, the degree of nonlinearity is 
typically moderate, i.e. both waves and eddies are present simultaneously. It is the relative 
importance of the anisotropic linear terms with respect to the isotropic nonlinear terms 
in the CHM equation which sets the anisotropy of the dynamics. If the linear terms are 
overpowered by the nonlinearity, the condensation of energy does not give rise to ZF's, but 
generate isotropic, round vortices. 

Related, but more simplified models beyond CHM, are the modified CHM^, Hasegawa- 
Wakatani (HW)^ model and modified Hasegawa-Wakatani modeUi. The HW model is given 
by equations ([1]) and ([2]) below. The term "modified" in reference to both the CHM and HW 
models means that the zonal-averaged component is subtracted from the electric potential 
to account for absence of the Boltzmann response mechanism for the mode which has no 
dependence in the direction parallel to the magnetic field. 

Quantitative investigations of the LH transition physics are presently carried out, using 
realistic modelling, such as gyrokinetic simulations, drawing inspiration from the qualitative 
results obtained by these idealised models. However, the understanding of the dynamics 
generated by these idealised models remains incomplete. It was only recently that the 
scenario of the drift- wave/ZF feedback loop proposed theoretically in 1990 for the CHM 
model was confirmed and validated by Direct Numerical Simulations (DNS) of the CHM 
equation by Connaughton et al^"^. In their work, the system was forced and damped by 
adding a linear term on the right-hand side of the CHM equation which would mimic a 
typical shape of a relevant plasma instability near the Larmor radius scales and dissipation 
at smaller scales. A drift-wave/ZF feedback loop was also seen in DNS of the modified HW 
model by Numata et al^^. The two dimensional simplification of the HW equations involves 
a coupling parameter generally called the adiabaticity. In one limit of this adiabaticity 
parameter the HW model becomes CHM model and in another limit it becomes the 2D 
Euler equation for an incompressible neutral fluid. The HW model contains more physics 
than CHM in that it contains turbulence forcing in the form of a (drift dissipative) instability 
and it predicts a non-zero turbulent transport - both effects are absent in CHM. On the other 
hand, it was claimed in Numata et al— that the original HW model (without modification) 
does not predict formation of ZF's. This claim appears to be at odds with the CHM results 
of Connaughton et al, considering the fact that HW model has CHM as a limiting case. 
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In the present work we will perform DNS of the HW model (without modification) aimed 
at checking realisability of the drift- wave/ZF feedback scenario proposed by Balk et al. 
in 1990^1^ and numerically observed by Connaughton et al — . This will be a step forward 
with respect to the CHM simulations because the instability forcing is naturally present in 
the HW model and there is no need to add it artificially as it was done for CHM. We will 
vary over a wide range of the coupling parameter of the HW model including its large values 
which bring HW close to the CHM limit. We will see that the ZF generation and turbulence 
suppression are indeed observed for such values of the coupling parameter, whereas for its 
smaller values these effects are lost. 



II. PHYSICAL MODEL 

The model we will consider is based on the HW equations^^: 

- X z ■ vj V^V" = e - n) - z/V^(VV), (1) 

- VV- X z ■ (n + ln{no)) = C - ra) - z/VV, (2) 

where ip is the electrostatic potential, n is the density fluctuation. The variables in Eq. ([1]) 
and Eq. ([2]) have been normalized as follows, 

x/ps ^ X , ujcit — )■ t , eip/Te — )■ , rii/riQ — )■ n. 



where ps = yTe/mu~^ is the ion gyroradius, no and ni are the mean and the fluctuating 
part of the density, e, m and Te are the electron charge, mass and temperature respectively 
and Uci is the ion cyclotron frequency. C is the adiabaticity parameter which we will discuss 
below. The last terms in Eq. ([1]) and Eq. ([2]) are 4*^-order hyperviscous terms that mimic 
small scale damping. 

The physical setting of the HW model may be considered as a simplification of the edge 
region of a tokamak plasma in the presence of a nonuniform background density uq = nQ{x) 
and in a constant equilibrium magnetic field B = -Bo^z, where is a unit vector in the 
z-direction. The assumption of cold ions and isothermal electrons allows one to find Ohm's 
law for the parallel electron motion: 

r]\\J\\ = -nev\i = E\\ + —V\\v = ^Vyn - Vy^, (3) 



where T^y is the parallel resistivity, E\\ is the parallel electric field, v\\ is the parallel electron 
velocity and p is the electron pressure. This relation gives the coupling of Eq. ([1]) and Eq. ([2]) 
through the adiabaticity operator C = —Tf./ {nQriujcie^)d'^ / dz^ . We will show in the following 
that the HW model describes the growth of small initial perturbations due to the linear 
drift-dissipative instability leading to drift wave turbulence evolving to generate ZF via an 
anisotropic inverse cascade mechanism followed by suppression of drift wave turbulence by 
ZF shear. 

The role of the dissipation terms in the equations ^ and ([2]) is to ensure the possibility 
of a steady state and to prevent a spurious accumulation of energy near the smallest resolved 
scales. We chose the dissipation terms proportional to uk*, but the qualitative picture is 
expected to be largely insensitive to the particular choice of the dissipation function. 

The important, relevant quantity in fusion research is the particle fiux in the x direction 
due to the fiuctuations. 



where k = ps|V/n(no)| is the normalized density gradient. Another quantity that we will 
monitor is the total energy, 



We will be interested in particular in the velocity field, the ZF's and their infiuence on 
turbulent fiuctuations. We therefore focus on the kinetic energy. Since one of the main 
subjects of the present work is the investigation of the ZF generation, we will quantify the 
energy contained in these ZF's by separating the energy into the kinetic energy Ez contained 
in a zonal sector, defined as \kx\ > \ky\, and the kinetic energy Em contained in a meridional 
SGctor, I I ^ I I * 

Different possibilities to determine the adiabaticity parameter will now be discussed. The 
first one is based on the linear stability analysis. In order to understand the dependence 
of the HW model on its constituent parameters it is useful to solve the linearised system. 
Linearisation of the Eq. ([1]) and Eq. ([2]) around the zero equilibrium {ip = and n = 0) and 
considering a plane wave solution, ip{k,t) ~ ■^^pe*'-'"'"'"'^*-* and n{k,t) ~ nQe^'^^'^^^^\ yields the 
dispersion relation for a resistive drift wave: 




(4) 




(5) 



cj^ + iu{h + 2vk'^) - ihu^ - Qvk'^{l + k"^) 



(6) 
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where = kyK/{l + is the drift frequency, b = C(l + k'^)/k'^ and k"^ = + ky. Let us 
introduce the real frequency ur and the growth rate 7 as 



(7) 



The dispersion relation ([6]) has two solutions, a stable one, with jmax = niax7(k) > 0, and 
an unstable one, with 7(k) < 0. The Fig. [TJ^a) shows the behaviour of 7 for k^ = and 
u = 10~^ for the unstable mode as a function of the wavenumber. The behaviour of 'jmax as 
a function of C is shown in Fig. [U^b). 
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FIG. 1. (a) Linear growth rate 7 for kx = 0, v = 10^^ and different values of C: C = 0.01 (solid 
line), C = 0.1 (dashed line), 6 = 1 (dotted line), C = 4 (dash-dot line), (b) Maximum linear growth 
rate jmax for fc^ = as a function of C, for u = 10~^. 



For the inviscid case, if u is ignored, the solution of Eq. 
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The maximum growth rate corresponds to 6 ~ 4a;^, and 
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(9) 



(1 + P)2' 

The adiabaticity operator C = —Te/{noricOcie'^)d'^/dz^ in Fourier space becomes an adia- 
baticity parameter via replacement d'^ /dz^ — > —k^, where k^ is a wavenumber characteristic 
of the fluctuations of the drift waves along the field lines in the toroidal direction. Recall 
that we assume the fluctuation length scale to satisfy the drift ordering, k\\ k±. It is 
natural to assume that the system selects kz which corresponds to the fastest growing wave 
mode. In this case one should choose the parallel wavenumber k^ which satisfies Eq. (Q 
for each fixed value of perpendicular wavenumber kj^ = {kx,ky). This approach is valid 



provided that the plasma remains colhsional for this value of kz- Note that according to 
Eq. IQ such a choice gives C = for the modes with ky = 0. 

Another common approach is to define the parameter C simply as a constant. This 
approach makes sense if the maximum growth rate correspond to values of kz which are 
smaller than the ones allowed by the finite system, i.e. kzmm = ^/R, where R is the bigger 
tokamak radius. In this case the HW model has two limits: adiabatic weak colhsional 
limit (C — i- oo) where the system reduces to the CHM equation, and the hydrodynamical 
limit (C — 7- 0) where the system of Eqs. ([T]) and reduces to the system of Navier-Stokes 
equations and an equation for a passive scalar mixing. 

In our simulation we will try and compare both approaches: choosing constant C and 
chosing 6 selected by the maximum growth condition (|9]). 

III. NUMERICAL METHOD 

Numerical simulations were performed using a pseudo-spectral Fourier code on a square 
box with periodic boundary conditions. The number of the modes varied from 256^ (with 
the lowest wavenumber Ak = 0.042 and the size of the box = Ly = 150) to 1024^ 
(with the lowest wavenumber Ak = 0.002 and the box size = Ly = 300), the viscosity 
coefficient was taken u = 0.0005 and u = 0.00005, respectively. The time integration was 
done by the fourth-order Runge-Kutta method. The integration time step was taken to be 
At = 5- IQ-^ and At = 10"^ 

IV. RESULTS 

In the following we will present our results for the evolution of the HW turbulence for 
different choices of the form and the values of the adiabaticity parameter. 

a. Constant adiabaticity In this section we will consider the case where the adiabaticity 
parameter is taken a constant. We present the results obtained with different values of C 
corresponding to the hydrodynamic regime (C — )■ 0, strongly colhsional limit), adiabatic 
regime (6 — )■ oo, weakly colhsional limit) and transition regime {C 1). The simulation 
parameters are presented in Table HI 

Fig. |2]shows the typical time evolution of the total kinetic energy, kinetic energy contained 
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TABLE I. The simulation parameters. 
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FIG. 2. Evolution of the total kinetic energy, energy contained in the zonal and the meridional 
sectors, (a) 6 = 0.01, (b) C = 1 and (c) C = 40. 

in the zonal sector \kx\ > \ky\, kinetic energy contained in the meridional sector \kx\ < \ky\ 
and the particle flux r„ (see Eq. (jlj)), for the adiabaticity parameter values C = 0.01, 1 and 
40, respectively. From Fig. [2] we see that the small initial perturbations grow in the initial 
phase. In this phase the amplitudes of the drift waves grow. Then, these drift waves start 
to interact nonlinearly. For the case C = 0.01 and 6 = 1 the resulting saturated state seems 
close to isotropic as far as can be judged from the close balance between and Em- For the 
simulation with C = 40 it is observed that the meridional energy strongly dominates until 
t ^ 4000. After this, the zonal energy rapidly increases and becomes dominant for t > 6000. 
This picture is in agreement with the scenario proposed in Connaughton et al^ for the 
CHM system. For the different values of C we can observe distinct types of behaviour in 
the evolution of the kinetic energy. The initial phase always agrees with the linear stability 
analysis (section [11]) • The speed at which the system enters to the saturated state is strongly 
dependent on C. The slowness of the transition of the system to a saturated level has limited 
the maximum value of the adiabaticity parameter to C = 40 and the maximum number of 
modes for such C to 256^. This slowness can be understood from the linear growth rate 
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dependence which decreases rapidly with C, see Fig. [U^b). For the value C = 0.01 and C = 1 
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(a) (b) 
FIG. 3. Fields of the stream functions, a) C = 1 t = 800, b) C = 40 t = 7600. 



we observe monotonous growth of the zonal, meridional and the total energies, as well as 
the the particle flux F^ - until these quantities reach saturation. We see that for C = 40 the 
initial growth of the meridional energy and the particle flux F„ is followed by a significant 
(between one and two orders of magnitude) suppression of their levels at the later stages. 
This is precisely the type of behavior previously observed in the CHM turbulence (Ref]l2l). 
and which corresponds to LH-type transport and drift-wave suppression. Recall that F„ is 
the particle flux in the x-direction which corresponds to the radial direction of the physical 
system that we model, the edge region of the tokamak plasma. 

Fig. [3](a) and Fig. [3](b) show instantaneous visualisations of the electrostatic potential ip 
for values 6 = 1 and C = 40, respectively. The structure of ip is strongly dependent on the 
regime: for low values of C the structure of is isotropic, whereas for the high values of 6 
the structure of ip is anisotropic and characterised by formation of large structures elongated 
in the zonal direction. 

For a better understanding of the anisotropic energy distributions, on Fig. H] we show the 
2D kinetic energy spectra normalized by their maxima for the cases 6 = 1 and C = 40. In 
the initial phase, for both cases, one can observe a concentration of the kinetic energy in the 
region corresponding to the characteristic scales of the drift wave instability, see Figs. Ht^a) 
and 111(b). Such a linear mechanism generates energy mainly in the meridional sector. For 
the saturated state, one can observe the distinct features of the energy distribution in the 
2D k-space. We can see that for large values of C there is a domination of concentration 
of the kinetic energy in the zonal sector, which absorbs energy from the meridional drift 



waves, see Fig. EJ^c). These computations for C ^ 1 are extremely long. Even though in 
the limit we should obtain the dynamics governed by the CHM equations^, this may only 
be approached for very high value C. Also the increase of C decreases the growth rate 
instability of the drift waves. Thus, comparison with Connaughton et aVs. simulation of 
CHM is not straightforward, since they artificially added a forcing term in order to mimic 
a HW-type instability and in the limit of C — )■ oo the HW system tends to the unforced 
CHM system. For C = 40 the zonal flows are not yet very pronounced in the physical space 
visuahzation Fig. Mjo), but very clear in Fourier space. Indeed, while in Fig. HJ^c) we see 
that the saturated 2D energy spectrum isotropic for C = 1, on Fig. IH^d) we can see that the 
spectrum is strongly anisotropic and mostly zonal for the C = 40 case. 
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FIG. 4. Snapshots of the kinetic energy spectrum normalized by its maximum value, (a) 6 = 1, 
transition regime; (b) C = 40, transition regime; (c) C = 1, saturated regime; (d) C = 40, saturated 
regime. 

b. Wavenumber dependent C. Now we will consider the case when the parameter C is 
defined according to the relation (Q for ps ~ 0.02 and k = 0.3491. For the given parameters 
maximum value of C equal to 0.453. Note that this case has in common with the MHW 
modelli that the coupling term in Eq. ([T]) and ([2]) is zero for the mode ky = 0. The 
numerical simulations were performed for z/ = 5 ■ 10^^ with 1024^ modes, box size L = 300 
and At = 10"^. Time evolution of the total, the zonal and the meridional kinetic energies, 
as well as the particle flux F^, are shown on Fig. [51 We observe a similar picture as before 



10 



in the simulation with large constant adiabaticity parameter C = 40. Namely, the total and 
the zonal energies grow monotonously until they reach saturation, whereas the meridional 
energy and the transport initially grow, reach maxima, and then get reduced so that their 
saturated levels are significantly less than their maximal values. This is because the ZF's 
draw energy from the drift waves, the same kind of LH-transition type process that we 
observed in the constant adiabaticity case with C = 40. In the final saturated state, there is 
a steady state of transfer of the energy from the drift waves to the ZF structures, so that the 
waves in the linear instability range in the meridional sector do not get a chance to grow, 
which can be interpreted as a nonlinear suppression of the drift-dissipative instability. 

10° 
10-1 
g 10-2 
; 10-3 

tij" 10-* 

10-5 
10-'' 

200 400 600 800 1000 
t 

FIG. 5. Evolution of the total energy (red line), energy contained in a zonal sector (green line), 
energy contained in a meridional sectors (blue line), particle flux F^ (magenta line). 

The physical space structure of the streamfunction ip is shown in Fig. |6]^a) and Fig. Mjo) 
for an early moment and for the saturated state. One can see formation of well-formed ZF's. 
Figs [71(a), (b) and (c) show the snapshots of the 2D energy spectrum evaluated at time 
t = 100, t = 190 and t = 800. We can see that initially meridional scales are excited via the 
linear instability mechanism, see Fig[71^a). This is followed by the nonlinear redistribution 
of the energy into the zonal sector, so that the spectrum for t = 190 looks almost isotropic, 
FigE(b)- The process of transfer to the zonal scales continues, and for t = 800 we observe 
a very anisotropic spectrum which is mostly zonal, see Fig[7](c). 

V. CONCLUSIONS 

In this paper we have studied numerically turbulence described by the Hasegawa- 
Wakatani model Eq. ([1]) and ([2]) for three different constant values of the adiabaticity 
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FIG. 6. Streamfunction field at (a) t = 190, (b) t = 800, for case C = e{k). 
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FIG. 7. Snapshots of the 2D kinetic energy spectrum normalized by its maximum value and 
evaluated at (a) t = 100, (b) t = 190 and (c) t = 800 for the case of k-dependent C. For figure (c) 
a close-up view of a zonal sector is shown in the right-hand corner. 

parameter C, and for a wavenumber-dependent C chosen to correspond to the fastest grow- 
ing modes of the drift-dissipative instabihty. Our aim was to resolve a visible contradiction 
between the assertion made by Numata et al}^ that zonal flows (ZF's) do not form in the 
original (unmodified) HW model and the clear observation of the ZF's by Connaughton et 
alr^ within the Charney-Hasegawa-Mima model which is a limiting case of the HW system 
for large 6. 

In our simulations for large values of 6, namely for C = 40, we do observe formation 
of a strongly anisotropic flow, dominated by kinetic energy in the zonal sector, followed by 
the suppression of the short drift waves, drift-dissipative instability and the particle flux, 
as originally proposed in the drift- wave/ZF feedback scenario put forward by Balk et al. 
in 1990^'^. This result suggests the original HW model to be the minimal nonlinear PDE 
model which can predict the LH transition. Note that even though the drift- wave/ZF loop 

12 



was also observed in the CHM simulations^^, it cannot be considered a minimal model for 
the LH transitions because the CHM model itself does not contain any instability, and it 
had to be mimicked by an additional forcing term. 

On the other hand, like in Numata et al.—, we see neither formation of ZF's nor sup- 
pression of the short drift waves and the transport for low values of C, 0.1 and 1. This is 
quite natural because in the limit of low C the HW model becomes a system similar to the 
isotropic 2D Navier-Stokes equation. Our guess is that Numata et al}^ have not explored 
the range of large C in their simulations, and therefore reached a conclusion that the origi- 
nal HW model would never produce ZF's. Note that the large C simulations are extremely 
demanding computationally because of the slow character of the linear instability in this 
case. 

For the simulation with wavenumber-dependent C, with a maximum parameter Gmax = 
0.453, we have also observed the predicted drift-wave/ZF feedback process characterised by 
even stronger and pronounced zonation than in the C = 40 case. Of course, a decision which 
case is more relevant, constant or wavenumber-dependent C, or the modified Hasegawa- 
Wakatani model^^ should be decided based on the plasma parameters and the physical 
dimensions of the fusion device. Namely, the wavenumber-dependent C should only be 
adopted if the fastest growing modes have wavenumbers allowed by the largest circumference 
of the tokamak; otherwise one should fix the parallel wavenumber at the lowest allowed value. 

In this paper once again we have demonstrated the usefulness of the basic PDE models of 
plasma, like CHM and HW, for predicting and describing important physical effects which 
are robust enough to show up in more realistic and less tractable plasma setups. Recall 
that the HW model is relevant for the tokamak edge plasma. As a step toward increased 
realism, one could consider a nonlinear three-field model for the nonlinear ion-temperature 
gradient (ITG) instability system, which is relevant for the core plasma; see eg. Leboeuf et 
aP^. This model is somewhat more complicated than what we have considered so far but 
still tractable by similar methods. This is an important subject for future research. 
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